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q Abstract 

X> 

We show that certain non-linear dynamical systems with non-linearities in the form of Hill 
i_i functions, can be approximated by piecewise linear dynamical systems. The resulting piecewise 

systems have closed form solutions that can be used to understand the behavior of the fully 
nonlinear system. We justify the reduction using geometric singular perturbation theory and 
^vq illustrate the results in networks modeling a genetic switch and a genetic oscillator. 

1 Introduction 

Accurately describing the behavior of interacting enzymes, proteins, and genes requires spatially 
extended stochastic models. However, such models are difficult to implement and fit to data, 

i-H hence modelers frequently use tractable reduced models. In most popular models of biological 

networks, the dynamics of each node is described by a single ODE, and sigmoidal functions are 

•'■h used to model interactions between the network elements. The resulting ODEs are generally 

not analytically tractable. This can hinder the study of large networks, where the number of 
parameters and the potential dynamical complexity make it difficult to analyze the behavior of 
the system using purely numerical methods. 

Analytical treatments are possible in certain limits. For instance, the approaches that have 
been developed to analyze models of gene interaction networks can be broadly classified into 
three categories [Polynikis et al. (2009)]: Quasi Steady State Approximations (QSSA), Piecewise 
Linear Approximations (PL A), and discretization of continuous time ODEs. 

Here, we aim to develop the theory of PL As. In certain limits interactions between network 
elements become switch-like [Kauffman (1969); Alon (2006); Davidich and Bornholdt (2008)]. 
For instance, the Hill function, f(x) = x n / (x n + J n ), approaches the Heaviside function, H (x — 
J), in the limit of large n. In this limit the domain on which the network is modeled is also 
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naturally broken into subdomains. For Hill functions, the thresholds, defined by J, divides the 
domain into two subdomains within which the Heaviside function is constant. Thus within each 
subdomain a node is either fully expressed, or not expressed at all. The original Hill function, 
f(x), is approximately constant in each of the subdomains, and boundary layers occur when x 
is close to threshold [Ironi et al. (2011)]. 

This general approach has a long and rich history, and piecewise linear functions of the form 
proposed in [Glass and Kauffman (1973)] have been shown to be well suited for the modeling of 
genetic regulatory networks (for a brief review see [De Jong (2002)]). In certain cases the results 
can be justified rigorously. In particular, singular perturbation theory can be used to obtain 
reduced equations within each subdomain and the boundary layers, and global approximations 
within the entire domain [Ironi et al. (2011)]. 

Here we take a similar approach, but work in a different limit. We again start with the 
Hill function, x n / (x n + J"), but assume that J is small. Although the subsequent results hold 
for any fixed n, for simplicity we assume n = 1. Equations involving this special class of Hill 
functions arc known as Michaelis-Menten equations, and J is known as the Michaclis-Menten 
constant [Michaelis and Menten (1913); Goldbeter and Koshland (1981); Ciliberto et al. (2007); 
Ma et al. (2009); Davidich and Bornholdt (2008); Goldbeter (1991); Novak and Tyson (1993); 
Novak et al. (2001); Tyson et al. (2003)]. We note that the models of chemical reactions we 
consider can be rigorously derived from the Chemical Master Equation only in the case of a 
single reaction [Kumar and Josic (2011)]. The models of networks of chemical reactions that we 
take as the starting point of our reduction should therefore be regarded as phenomenological. 

We will examine the case when the Michaelis-Menten constant, J, is small. This case has a 
simple physical interpretation: Consider the Hill function that occurs in the Michaelis-Menten 
scheme, where an enzyme is catalyzing the conversion of the inactive form of some protein to 
its active form. When J is small the total enzyme concentration is much smaller than the 
total protein concentration. The asymptotic limit J — > was recently considered to obtain 
heuristically a Boolean approximation of a protein interaction network [Davidich and Bornholdt 
(2008)]. Here we consider a rigorous justification underlying such reductions, as well as how the 
reduction could be used to understand the dynamics of gene networks. 

The main idea behind the reduction we propose can be summarized as follows: Given the 
non-linear term f(x) = x/(x + J), when x >• J then f(x) ~ 1, and when x m then we do the 
analysis by introducing a new variable like x :— J/x. This new variable x serves as a microscope 
to observe the boundary regions. As we will show, the domain is naturally decomposed into 
a nested sequence of hypercubes such that for each level of nesting we get a separate linear 
equation. 

We proceed as follows: In Section 2 we illustrate our approach using simple examples and 
provide numerical evidence for the validity of our claim. In Section 3 we describe a general class 
of differential equations which subsumes these examples. Furthermore, in this section we justify 
our approach mathematically using Geometric Singular Perturbation Theory (GSPT). We will 
conclude with a discussion on limitations of these reductions. 

2 Example problems 

We start by demonstrating the main idea of our approach in the cases of two and three mutually 
repressing biological elements. For instance, these elements could be genes that mutually inhibit 
each other's production [Gardner et al. (2000); Elowitz and Leibler (2000)]. However, as the 
theory we develop is general, we do not constrain it to a particular interpretation. We first 
provide an intuitive illustration of the approach along with a heuristic justification of the different 
steps in the reduction. A mathematical justification follows. 
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(a) toggle switch (b) repressilator 

Figure 1: (a) Nodes ui, U2 inhibiting each others activity. End result is like a switch. The node which 
was stronger in the beginning will stay stronger and will completely suppress the other, (b) Nodes Ui, U2, 
and it3 suppressing each other in a cyclic fashion. Not surprisingly, the end result is oscillatory behavior. 



2.1 A network of two mutually inhibiting elements 

First we consider two mutually repressing elements within a biological network. This toggle 
switch motif (see Figure la) is common in biological networks [Tyson et al. (2003); Gardner 
et al. (2000)]. Let ui,u 2 G [0,1] represent the normalized levels of activity of the first and 
second element, respectively. Therefore, when Ui — 1 the i th network element is maximally 
active (expressed). The system be modeled by 

dux „ _ 1 - Ui ui 

—— = 0.5 U2 , 

at J + 1 — ui J + ui , . 

du 2 1 - u 2 u 2 

— — = 0.5 ui , 

at J + 1 — u 2 J + u 2 

where J is some positive constant. The structure of Eq. (1) implies that the cube [0, l] 2 = 
{(ui,u 2 ) | < Ui,u 2 < 1} is invariant (see Proposition 1). 

In the limit of small J, Eq. (1) can be approximated by a piecewise linear differential equation 
as follows: If itj is not too close to zero the expression Ui/{J + ut) is approximately unity. More 
precisely, we fix a small 6 > 0, which will be chosen to depend on J. When m > S and J is 
small then m/ (J + Ui) ~ 1. Similarly, when ui > 1 — 5 then (1 — Uj)/(J + 1 — Ui) 1. 

With this convention in mind we break the cube [0, l] 2 into several subdomains, and define 
a different reduction of Eq. (1) within each. For example, the interior of the domain [0, l] 2 is 
defined by 

Kl := {(ui, u 2 ) e [0, l] 2 | 5 < ui < 1 - 6 and S < u 2 < 1 - 6}. (2) 
Eq. (1), restricted to TZ® is approximated by the linear differential equation 

dui nK du 2 

_ =0 .5- W2 , _=0.5- U i. (3) 

On the other hand, if one of the coordinate is near the boundary, while the other is in the 
interior, the approximation is different. For instance, the region 

Kl := {{ui,u 2 ) e [0,1] 2 \ui < S and S< u 2 < 1 -6}, (4) 

forms a boundary layer where u\ is of the same order as J. Therefore the term Ui/{J + Ui) can 
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not be approximated by unity. Instead the approximation takes the form 



— = 0.5-u 2 — , (5a 

at J + ui 

^=0.5-„„ (5b) 

This equation can be simplified further. Since the boundary defined by u\ — is invariant, 
must be small inside the boundary layer 72.J. We therefore use the approximations ^ w in 
Eq. (5a) and Ui w in Eq. (5b) to obtain 

= 0.5-u 2 Ml , (6a) 

J + Ul 

£-a* m 

Note that Eq. (6b) is linear and decoupled from Eq. (6a), while Eq. (6a) is an algebraic system 
which can be solved to obtain u\ w J/(2u2 — 1). Within TZ® we thus obtain the approximation 
u 2 (t) « 0.5t + u 2 (0) and m(t) w J/(t + 2u 2 (0) - 1). 

Note that here we have the freedom of only specifying the initial condition u 2 (0), while 
zti(O) is determined from the solution of the algebraic equation (6a). As we explain below, this 
algebraic equation defines a slow manifold within the subdomain The reduction assumes 
that solutions arc instantaneously attracted to this manifold. 

Table 1 shows how these ideas can be extended to all of [0, l] 2 . In each of the 9 listed 
subdomain one or both variables are close to either or 1. Therefore each subdomain corresponds 
to either the interior, edge, or corner of the unit square. Following the preceding arguments, 
we assume that variable(s) that are close to or 1 are in steady state and lead to an algebraic 
equation. Similarly, the evolution of the interior variables is described by linear differential 
equations. The resulting algebraic-differential systems are given in the last column of Table 1. 

The reductions in the corner subdomains 1Z® 2 , 7\Ljy 2 , 1Z\, and 1Z\ consist of purely algebraic 
equations. When J is small some of these equations will have a solution in [0, l] 2 , indicating a 
stable fixed point near the corresponding corner (Jl\ and 1Z\)- Others will not have a solution in 
[0, l] 2 , indicating that approximate solutions do not enter the corresponding subdomain (R\_ 2 
and 7\Lq :2 ). 

Each approximate solution has the potential of exiting the subdomain within which it is 
defined, and entering another. The global approximate solution of Eq. (1) is obtained by using 
the exit point from one subdomain as the initial condition for the approximation in the next. In 
subdomains other than TZq some of the initial conditions will be prescribed by the algebraic part 
of the reduced system. The global approximation may therefore be discontinuous, as solutions 
entering a new subdomain are assumed to instantaneously jump to the slow manifold defined 
by the algebraic part of the reduced system. Fig. 2 shows that when J is small, this approach 
provides a good approximation. 

2.2 A network of three mutually inhibiting elements 

The same reduction can be applied to systems of arbitrary dimension. As an example consider 
the repressilator [Tyson ct al. (2003); Elowitz and Leiblcr (2000)] described by 

dui 1 - ui ui 
0.6— — u 3 - 



dt J + 1 — ui J + ui ' 

du 2 1 - u 2 u 2 

—r~ = 0.4— ui— , (7) 

dt J + 1 - u 2 J + u 2 w 

du 3 1 - u 3 u 3 

—r- = 0.3— U 2 - . 

at J + 1 — u 3 J + u 3 
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Figure 2: Comparison of the numerical solution of Eq. (1) (dashed black) and the solution of the 
approximate system as listed in Table 1 (solid colored) for two different values of J ( We used J = 10~ 2 
in (a); and J = 10~ 4 in (b).). The different colors denote the switching behavior of the solution from one 
subdomain to next. We used S = 0.01. Solution of the linear approximation started in the subdomain 
TZq (Initial value: ui — 0.6, U2 = 0.4), and as soon as U2 became smaller than 5, the subdomain switched 
to TZ% and driving linear differential equation also switched accordingly. It should be noted that the 
approximate solution is discontinuous. The reason is that as soon as the solution crossed the horizontal 
line, «2 = S, the solution jumped (see inset) to the manifold, described by the algebraic part of the 
linear differential algebraic system prevalent in the subdomain 1Z®. The solution finally stopped in the 
subdomain IZ2. 



The cyclic repression of the three elements in this network leads to oscillatory solutions over a 
large range of values of J. The domain of this system, [0, l] 3 , can be divided into 27 subdo- 
mains: 1 interior, 6 faces, 12 edges, and 8 vertices. We can again approximate Eq. (7) with 
solvable differential-algebraic equation within each subdomain, to obtain a global approximate 
solution. We demonstrate the validity of this approximation in Fig. 2.2. Note that both the 
numerically obtained solution to Eq. (7), and its approximation exhibit oscillations, and that 
the approximation is discontinuous. 



3 General setup 

The approximations described in the previous section can be extended to more general models. 
Suppose we describe the evolution of n interacting elements, u\, U2, u n , by 



dui , 1 — Ui 



-Wt^T. (8) 



dt 1 Jf + l-Ui 1 Jf + u. 

where J/ 1 , J/ are some positive constants. Here A$ and Ii are activation/inhibition functions 
that capture the impact of other variables on the evolution of itj The initial conditions are 
assumed to satisfy Uj(0) € [0, 1] for all i. 

We assume that the activation and inhibition functions are both affine [De Jong (2002)], 

n n 

j=i 3=1 

where we use the convention x + — max{a;,0} and x~ = max{— x, 0}. The n x n matrix, 
W — [wij] and the n x 1 vector b = [b\ hi ... b n ] capture the connectivity and external input 
to the network, respectively. In particular, gives the contribution of the j th variable to the 
growth rate of i th variable. If Wij > 0, then uiy appears in the activation function for Ui\ and 
if Wij < then — w^j appears in the inhibition function for itj. The intensity of the external 
input to the i th element is \bi\, and it contributes to the activation or the inhibition function, 
depending on whether bi > or hi < 0, respectively. 
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Figure 3: Comparison of the numerical solution of Eq. (7) (dashed black) and the solution of the 
approximate linear system (not explicitly provided) for two different sets of J and S. For (a)-(c) J = 
10 _2 ,<5 = 0.06; for (c)-(f) J = 10 _4 ,<5 = 0.01. The approximate solution changes color when switching 
between different subdomains. Note that the approximate solution is discontinuous in general. The 
reason is that as soon as the solution enters a new subdomain, the solution jumps (see inset) to the 
manifold defined by the algebraic part of the linear differential algebraic system corresponding to the new 
subdomain. 



Proposition 1. If Ai and Ii are positive, then the cube [0,1]™ is invariant for the dynamical 
system given by Eq. (8). 

Proof. It will be enough to show that the vector field at any point on the boundary is directed 
inward. Since, Ai and Ii are positive, for any i, 



dui 
~dt 



= A. 



1 



u i= 



1 tA 



> 0, and 



duj 
dt 



= -I 



1 



< o. 



□ 



4 General reduction of the model system 

To obtain a solvable reduction of Eq. (8) we follow the procedure outlined in Section 2. We 
present the result here, and provide the mathematical justification in the next section. For 
notational convenience we consider the case = J- = J, with J small and positive. The 
general case is equivalent. Let S be some positive number which will be used to define the 
thickness of the boundary layers, and which will depend on J in general. We start with the 
subdivision of the n-dimensional cube, [0, 1]™. 

Let T and S be two disjoint subsets of {1, 2, n}, and let 



K T S 



Uu u u 2 ,...,u„) G [0,1]™ 



u s < S for all s G S; u t > 1 — 5 for all t G T; 

and 5 <u k < 1 - S for all fc^Surj. 
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We extend the convention used in Table 1, and in Eqs. (2) and (4) so that TZq := lig when S 
is empty; IZg :— IZg when T is empty; and TZq := T&S wnen T, S are both empty. 

Within each subdomain TZ^ Eq. (8) can be approximated by a different linear differential- 
algebraic system. Following the reduction from Eq. (1) to Eq. (5), for i £ S U T we obtain the 
linear system 

dui \ 

-fa=Z*, a H u i+ b i- ( 10a ) 
J=l 

For s € 5 one of the nonlinear terms remains and we obtain 

f = (±<u, + tf) - (|>^ + 7T^' < 10b > 

while for i e T we will have 

£ = (g^ + tf) - + >r) ■ (loo 

Eq. (10) is simpler than Eq. (8), but it is not solvable yet. Following the reduction from Eq. (5) 
to Eq. (6), we now further reduce Eqs.(lOb-lOc). First we use the approximations u s rts and 
u t ~ 1 in the activation and inhibition functions appearing in Eq. (10). Second, we assume that 
u s for s e S and u t for t € T are in steady state. 

Under these assumptions we obtain the reduction of Eq. (8) within any subdomain TZ^ 

^ = E a ^ u 3 + ^2 a ij+ b -' i ^ S UT; 

j£SUT jer 

j£SUT t£T \jgSUT t£T J s 

(lib) 

° = - E ^i+E^-H' + E w+E a y +b *"' <eT - 

(lie) 

Eq. (11) is solvable since Eq. (11a) is decoupled from the rest, and Eqs. (lib) and (11c) are 
solvable for u s and u t , respectively, as functions of the solution of Eq. (11a). 

5 Mathematical justification 

We next justify our claim that the variables close to the boundary can be assumed to be in 
steady state. We define the following new variables to "magnify" the boundary region. 

m s := y for s e S, and u t := ^—p- for t G T. (12) 
Using Eq. (12) in Eq. (10) we get for i <£ S U T 



dui 
~dt 



dijUj + E + ^ I E ais " s " E ai *"* J + bu ( 13a ) 
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and for s e S, 

jdu^_ = j2 a + jUj +J2 a tt + J \ J2 a % ^ ~ E a ^" 4 + b > 
jgsur teT \jes teT J 

E w + E a ^ + - ■/ E a i-^ - E a ^* 5 ttV' (13b) 

V^SUT teT / i + Us \ieS teT / 1 + 

and similarly, for t € T, 



J ^r =-[ E a *>j ■ + E4 + b t ytV " J E a * + ^ -E4 fi i r 



+ E + E a *j ■ + b t + J E a i" s ~ E a «"j • ( 13c ) 

jgSUT jET \seS iST / 

When J is small, we can apply Geometric Singular Perturbation Theory (GSPT) to Eq. (13) 
[Hek (2010); Kaper (1998)]. The GSPT posits that, under a normal hyperbolicity condition 
which we will prove below, Eq. (13) can be further simplified by assuming that J = 0. This 
yields a differential-algebraic system 

= E a-ijUj + y^a>ij +bj, i^SUT: 
itSUT jet 

(14a) 

°= E a ij u o+J2 a tt+ b t- [ E a si u J +J2 a *t+K ) r^V' se5 ; 

jgSUT teT \i^SuT teT y 8 

(14b) 



° = -( E_4^+E<+ 6 * + ) E_w + E a «+ 6 *"' teT - 

(14c) 



\j<£SUT ]£T / ' Ut jgSUT ]ET 



which is equivalent to Eq. (11) after rescaling. This conclusion will be justified if the manifold 
defined by Eqs. (14b) and (14c) is normally hyperbolic and stable [Fenichel (1979); Kaper (1998); 
Hek (2010)]. We verify this condition next. 

Let u = {u il , ...,Ui m } where ...,i m } = {1,2,..., n}\(SUT), be the coordinates of u which 
are away from the boundary, and denote the right hand side of Eq. (14b) by F s (u,u is ), for all 
s e S, so that 

F s (u,u ie ):= a tj u 3+^2 a tt + b t ~ E a sj u 3 +^2 a st + K J 

j(SuT teT \i^SuT teT J 

and 

2 



+ U s 



= ( E a; jUj +J2a- t + b^j (y^t) ><° 



du 



for all s e S. Similarly, by denoting the right hand side of Eq. (14c) by G t (u,u it ), for all t e T. 



i.e. 



G t (u, u it ) := - I a tj u j + E a tj + b t Y+u~ + E a tj u i + E a tj + & * ' 

U^sut jeT y * j^suT jeT 
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we see that 



dG t 



( 



i6T 



) 




l + u t 



) 



2 



< 0. 



Hence, the manifold defined by Eqs. (14b) and (14c) is normally hyperbolic and stable. This 
completes the proof that the reduction of the non-linear system (8) to a solvable system (11) is 
justified for small J. 

6 Discussion 

A special class of non-linear differential equation was studied with non-linear interaction terms 
given by Hill functions. We showed that when the Michaelis-Mcnten constants are sufficiently 
small, the behavior of the system is captured by an approximate piecewise linear systems. 
This induces a natural decomposition of the domain into a nested sequence of hypercubes, 
with a separate linear-algebraic system giving an approximation in each subdomain. We have 
illustrated the theory in examples, and justified the conclusions using GSPT. 

A potential limitation in our arguments is that we have an approximation valid only in an 
asymptotic limit. It is unknown when and how the approximation breaks down. Another major 
limitation of our analysis is that we have not provided a systematic relationship between the 
thickness of the boundary, 5, and the Michaelis-Menten constant, J. Numerical tests suggest 
that J = 0(5 2 ). 

References 

Uri Alon. An Introduction to Systems Biology: Design Principles of Biological Circuits. Chap- 
man and Hall/CRC, 1 edition, July 2006. 

Andrea Cilibcrto, Fabrizio Capuani, and John J. Tyson. Modeling networks of coupled en- 
zymatic reactions using the total quasi-steady state approximation. PLoS Computational 
Biology, 3(3)(3):e45, March 2007. 

M. Davidich and S. Bornholdt. The transition from differential equations to boolean networks: 
A case study in simplifying a regulatory network model. Journal of Theoretical Biology, 255 
(3):269-277, December 2008. 

H. De Jong. Modeling and simulation of genetic regulatory systems: a literature review. Journal 
of computational biology, 9(1):67-103, 2002. 

Michael B. Elowitz and Stanislas Leibler. A synthetic oscillatory network of transcriptional 
regulators. Nature, 403(6767):335-338, January 2000. 

N. Fenichel. Geometric singular perturbation theory for ordinary differential equations. Journal 
of Differential Equations, 31(l):53-98, January 1979. 

Timothy S. Gardner, Charles R. Cantor, and James J. Collins. Construction of a genetic toggle 
switch in escherichia coli. Nature, 403(6767) :339-342, January 2000. 

L. Glass and S. Kauffman. The logical analysis of continuous, non-linear biochemical control 
networks. Journal of Theoretical Biology, 39(1):103-129, April 1973. 

A. Goldbeter. A minimal cascade model for the mitotic oscillator involving cyclin and cdc2 
kinase. Proceedings of the National Academy of Sciences of the United States of America, 88 



(20):9107-9111, October 1991. 



9 



Albert Goldbeter and Daniel E. Koshland. An amplified sensitivity arising from covalent modi- 
fication in biological systems. Proceedings of the National Academy of Sciences of the United 
States of America, 78(ll):6840-6844, 1981. 

Geertje Hek. Geometric singular perturbation theory in biological practice. Journal of mathe- 
matical biology, 60(3):347-386, March 2010. 

Liliana Ironi, Luigi Panzeri, Erik Plahte, and Valeria Simoncini. Dynamics of actively regulated 
gene networks. Physica D, 240:779-794, 2011. 

Tasso J. Kaper. An introduction to geometrical methods and dynamical dystems for singu- 
lar perturbation problems. In Analyzing Multiscale Phenomena Using Singular Perturbation 
Methods: American Mathematical Society Short Course, January 5-6, 1998, Baltimore, Mary- 
land (Proc. Sym. Ap.), pages 85-132, 1998. 

S.A. Kauffman. Metabolic stability and epigenesis in randomly constructed genetic nets. Journal 
of theoretical biology, 22(3):437-467, 1969. 

Ajit Kumar and Kresimir Josic. Reduced models of networks of coupled enzymatic reactions. 
Journal of Theoretical Biology, March 2011. 

Wenzhe Ma, Ala Trusina, Hana El-Samad, Wendell A. Lim, and Chao Tang. Defining network 
topologies that can achieve biochemical adaptation. Cell, 138(4) :760-773, August 2009. ISSN 
00928674. doi: 10.1016/j.cell.2009.06.013. 

Leonor Michaelis and Maud Menten. Die kinctik der inwertin wirkung. Biochemische Zeitschrift, 
(49):333-369, 1913. 

B. Novak and J. J. Tyson. Numerical analysis of a comprehensive model of m-phase control 
in xenopus oocyte extracts and intact embryos. Journal of cell science, 106(4)(4):1153-1168, 
December 1993. 

Bela Novak, Zsuzsa Pataki, Andrea Ciliberto, and John J. Tyson. Mathematical model of the 
cell division cycle of fission yeast. Chaos (Woodbury, N. Y.), ll(l):277-286, March 2001. 

A. Polynikis, S. J. Hogan, and M. di Bernardo. Comparing different ODE modelling approaches 
for gene regulatory networks. Journal of theoretical biology, 261(4):511-530, December 2009. 

J. J. Tyson, K. C. Chen, and B. Novak. Sniffers, buzzers, toggles and blinkers: dynamics of 
regulatory and signaling pathways in the cell. Current Opinion in Cell Biology, 15(2):221-231, 
2003. 



10 



Subdomain's name 




U2 


Approximating linear system 


K 


6<U!<1-S 


5 < u 2 < 1 - 5 


u[ = 0.5 — U2, 
u 2 = 0.5 — ui 


nl 


ui > 1 - 6 


5<u 2 <l-5 


1 — Ui 

= 0.5 -u 2 , 
J + 1 — U\ 

u' 2 = -0.5 


nl 


5 < ui < 1-5 


U2 > 1 — 5 


u[ = -0.5, 

1 — Mo 

= 0.5 T iti 

J + 1 - U 2 




U\ < 


< U2 < 1 — 


= 0.5 -u 2 — 

J + Ul 

v! 2 = 0.5 


n\ 


5 < ui < 1-5 


U 2 < 5 


u[ = 0.5, 
= 05 ui U2 

J + U 2 


nf 


ui > 1 - 5 


U2 > 1 — 5 


1 — Ul 

= 0-5 1 -1, 

J + 1 — Ui 

n - 1_U2 1 

U — U.O 1 

J + 1 - ita 


n\ 2 


ui < 5 


U2 < 5 


= 0.5 J T U1 , 

J + Ui 

= 0.5 J/ 2 

J + u 2 


n\ 


ui > 1 - 6 


U2 < 5 


1 — Ul 

= 0.5 

J + 1 — Ul 

= 0.5 / 2 

J + u 2 


n\ 


u\ < 5 


u 2 > 1 - 5 


= 0.5- ^ , 

J + Ul 

1 — U 2 
= 0.5- — 

J + 1-U2 



Table 1: List of differential-algebraic systems that approximate Eq. (1) in different parts 
of the domain. The subdomains are named so that the superscript (subscript) lists the co- 
ordinates that are close to f (close to 0), with denoting the empty set. For example, 1Z\ 
denotes that subdomain with u\ « 1 and ui ^ 0, and TZq the subdomain where ui is near 
1, but Mi is away from the boundary. The middle column define the subdomain explicitly. 
The right column gives the differential-algebraic system that approximates Eq. (1) within the 
given subdomain. 
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